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Many  processes  important  to  chemistry,  materials  science,  and  biology  cannot  be  described  without 
considering  electronic  and  nuclear-level  dynamics  and  their  coupling  to  slower,  cooperative  motions 
of  the  system.  These  inherently  multiscale  problems  require  computationally  efficient  and  accurate 
methods  to  converge  statistical  properties.  In  this  paper,  a  method  is  presented  that  uses  data  directly 
from  condensed  phase  ab  initio  simulations  to  develop  reactive  molecular  dynamics  models  that 
do  not  require  predefined  empirical  functions.  Instead,  the  interactions  used  in  the  reactive  model 
are  expressed  as  linear  combinations  of  interpolating  functions  that  are  optimized  by  using  a  lin¬ 
ear  least-squares  algorithm.  One  notable  benefit  of  the  procedure  outlined  here  is  the  capability  to 
minimize  the  number  of  parameters  requiring  nonlinear  optimization.  The  method  presented  can  be 
generally  applied  to  multiscale  problems  and  is  demonstrated  by  generating  reactive  models  for  the 
hydrated  excess  proton  and  hydroxide  ion  based  directly  on  condensed  phase  ab  initio  molecular  dy¬ 
namics  simulations.  The  resulting  models  faithfully  reproduce  the  water-ion  structural  properties  and 
diffusion  constants  from  the  ab  initio  simulations.  Additionally,  the  free  energy  profiles  for  proton 
transfer,  which  is  sensitive  to  the  structural  diffusion  of  both  ions  in  water,  are  reproduced.  The  high 
fidelity  of  these  models  to  ab  initio  simulations  will  permit  accurate  modeling  of  general  chemical  re¬ 
actions  in  condensed  phase  systems  with  computational  efficiency  orders  of  magnitudes  greater  than 
currently  possible  with  ab  initio  simulation  methods,  thus  facilitating  a  proper  statistical  sampling 
of  the  coupling  to  slow,  large-scale  motions  of  the  system.  ©  2012  American  Institute  of  Physics. 
[http://dx.doi.org/10.1063/L4743958] 


I.  INTRODUCTION 

Many  phenomena  in  condensed  phase  systems  of  interest 
to  chemists,  biologists,  and  material  scientists  involve  chem¬ 
ical  processes  spanning  multiple  time  and  length  scales.  This 
situation  introduces  various  challenges  that  need  to  be  over¬ 
come  if  one  is  to  simulate  these  systems  with  a  sufficiently 
accurate  model  to  capture  the  essential  physics  and  prop¬ 
erly  sample  the  required  system  sizes  and  time  scales  in  or¬ 
der  to  statistically  converge  the  calculated  properties  of  in¬ 
terest.  Condensed  phase  simulation  methods  that  explicitly 
treat  electronic  degrees  of  freedom, which  include  the  pos¬ 
sibility  for  modeling  chemical  reactions,  are  naturally  the  first 
method  of  choice;  but  they  come  at  a  significant  computa¬ 
tional  expense,  thus  limiting  the  accessible  length  and  time 
scales.  Also,  current  ab  initio  molecular  dynamics  (AIMD) 
simulations  that  are  calculated  by  using  density  functional 
theory  require  validation  from  more  accurate  Schrodinger 
wavefunction  methods.  Empirical  models  that  retain  the  accu¬ 
racy  and  essential  physics  of  electronic  structure  methods  can 
instead  be  used  to  calculate  molecular  simulations  at  a  greatly 
reduced  computational  cost,  thus  making  them  invaluable  for 
modeling  phenomena  in  heterogeneous  condensed  phase  sys¬ 
tems,  which  may  contain  several  long  length  and  timescale 
components  and  complex  interfaces.  The  situation,  however. 
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requires  more  care  when  chemical  reactions  are  involved; 
and  an  accurate  determination  of  the  reaction  barriers  in  the 
condensed  phase  environment  is  necessary.  The  latter  also 
requires  a  proper  statistical  sampling  of  slow  degrees  of  free¬ 
dom  that  may  be  coupled  to  the  reaction  of  interest.  The  mul¬ 
tistate  empirical  valence  bond  (MS-EVB)  method'1-7  is  one 
such  empirical  reactive  MD  method  that  has  been  exclusively 
used  to  model  the  complex  process  of  proton  solvation  and 
transport. s 

The  accurate  modeling  of  a  condensed  phase  system 
is  particularly  challenging  when  multiple  electronic  excited 
states  are  involved  as  typical  of  nonadiabatic  processes.  The 
computational  cost  for  the  nonadiabatic  ab  initio  simulation 
methods  can  be  more  expensive  than  ground-state  methods, 
which  further  limits  the  applicable  system  sizes  and  times 
that  can  be  simulated.912  Highly  accurate  methods  such  as 
coupled-cluster  algorithms13  are  sometimes  required  to  prop¬ 
erly  describe  excited  states,  but  these  methods  are  currently 
limited  to  small  molecules,  thus  making  any  attempt  at  sim¬ 
ulating  nonadiabatic  chemical  reactions  in  condensed  phase 
systems  with  the  same  method  infeasible.  Clearly,  a  suffi¬ 
ciently  general  algorithm  is  needed  to  develop  reactive  MD 
models  (adiabatic  or  nonadiabatic)  that  are  independent  of 
the  electronic  structure  method  used  for  the  reference  calcu¬ 
lations  and  that  are  able  to  account  for  variable  bond  topol¬ 
ogy  to  model  chemical  reactions.  These  models  could  then  be 
used  to  address  the  challenges  of  converging  statistical  prop¬ 
erties  relating  to  multiscale  reactive  phenomena  in  complex, 
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condensed  phase  systems.  The  work  described  in  this  paper 
provides  such  a  model. 

Several  MS-EVB  models  have  been  recently  developed 
to  model  protonic  defects. 14-19  Charge-ring  models  have  suc¬ 
cessfully  reproduced  many  structural  properties  of  aqueous 
hydroxide,16,20,21  but  MS-EVB  models  with  this  architecture 
have  been  shown  to  suffer  from  large,  currently  undiagnosed, 
energy  drifts.16  Wick  and  Dang13  and  Wick18  used  experi¬ 
mental  properties  and  ab  initio  calculations  of  clusters  to  pa¬ 
rameterize  MS-EVB  models  for  hydroxide  and  hydronium. 
The  hydronium  model18  results  in  a  more  localized  excess 
proton  (ci2  =  0.896  ±  0.002)  when  compared  to  the  previ¬ 
ously  reported  occupation  of  MS-EVB  states  (ci2  ~  0.6). 19 
The  increased  localization  of  the  excess  proton  will  result  in 
a  reduced  effective  radius,  which  is  expected  to  affect  defect 
transport  and  affinity  for  the  liquid-vapor  interface.  Some  of 
these  models,15,18  including  one  for  bulk  water,22  incorporate 
molecular  polarization  via  configuration  dependent  electro¬ 
statics  (e.g.,  bond  distance-dependent  atomic  charges  and/or 
induced  point  polarizabilities),  in  addition  to  the  delocaliza¬ 
tion  of  the  proton  defect  across  multiple  molecules  in  MS- 
EVB,  which  may  lead  to  an  unphysical  overestimation  of  po¬ 
larization.  To  our  knowledge  there  have  been  no  extensive 
studies  of  the  charge  distribution  of  molecular  models  con¬ 
taining  both  configuration  dependent  electrostatics  in  multi¬ 
state  reactive  models,  but  AIMD  simulations  and  nonreactive 
polarizable  simulations  have  been  shown  to  result  in  distinctly 
different  induced  dipole  distributions.23  Instead  of  parameter¬ 
izing  models  based  on  gas-phase  clusters,  this  work  aims  to 
exclusively  employ  information  from  condensed  phase  ab  ini¬ 
tio  simulations  to  systematically  force  match  reactive  models 
based  on  results  in  a  self-consistent  forcefield  faithful  to  the 
reference  ab  initio  simulation. 

In  our  previous  work,  a  force-matching  (FM)  algorithm 
was  developed  that  used  data  directly  from  a  condensed  phase 
ab  initio  simulation  to  determine  effective  interactions  be¬ 
tween  nuclei  for  an  empirical  reactive  MS-EVB  model  of  the 
hydrated  excess  proton.24  The  FM  method  was  sufficiently 
flexible  that  those  interactions  expected  to  dominate  far  from 
the  reactive  species  could  be  defined  by  using  previously  pa¬ 
rameterized  empirical  models  for  those  regions  (e.g.,  pure  sol¬ 
vent  or  protein  forcefield).  Those  portions  of  the  model  de¬ 
scribing  interactions  with  and  between  the  reactive  species 
could  then  be  separately  fit  by  using  data  from  AIMD  simu¬ 
lations.  However,  in  that  earlier  work  the  functional  forms  for 
the  interactions  in  the  reactive  MD  model  were  chosen  to  be 
the  same  predefined  empirical  functions  that  had  been  used  in 
an  earlier  MS-EVB  model  for  the  excess  proton.19  The  algo¬ 
rithm  presented  here,  on  the  other  hand,  provides  a  significant 
advantage  by  enabling  the  linear  least-squares  optimization 
of  tabulated  expressions,  thereby  removing  any  possible  ar¬ 
tifacts  introduced  by  empirical  functions  that  may  constrain 
the  accuracy  of  the  resulting  FM  model.  Linear  FM  methods 
also  have  the  advantage  of  no  longer  requiring  the  nonlinear 
optimizations  of  several  parameters,  which  may  be  slow  to 
converge  or  plagued  with  multiple,  possibly  unphysical,  lo¬ 
cal  minima.  The  present  model  can  therefore  be  considered  to 
be  “multiscale”  in  that  it  connects  electronic  structure  to  ef¬ 
fective  reactive  forces  between  the  nuclei,  but  it  is  no  longer 


“empirical”  in  the  sense  that  predefined  empirical  analytical 
functions  are  not  used  to  represent  those  effective  forces. 

The  procedure  used  to  develop  these  new  reactive  MD 
models  is  actually  a  generalization  of  the  multiscale  coarse- 
graining  (MS-CG)  method.25-29  The  MS-CG  method  reduces 
the  complexity  of  the  original  system  by  constructing  a  re¬ 
duced  coarse-grained  system  with  significantly  fewer  degrees 
of  freedom,  the  dynamics  of  which  are  governed  by  a  set 
of  effective  interactions.  For  AIMD  simulations,  the  original 
nuclear+electron  system  can  thus  be  mapped  to  a  reduced 
representation  containing  only  nuclei.21  With  the  electronic 
degrees  of  freedom  integrated  away,  a  set  of  effective  in¬ 
teractions  between  nuclei,  determined  by  using  a  variational 
force-matching  algorithm,  can  be  used  to  sample  the  long¬ 
time  motion  of  the  nuclear  degrees  of  freedom  at  a  signif¬ 
icantly  reduced  computational  cost.  The  elimination  of  the 
electronic  degrees  of  freedom  is  analogous  to  the  develop¬ 
ment  of  “solvent-free”  coarse-grained  models  in  which  the 
solvent  degrees  of  freedom  are  integrated  away.30-'1  The  re¬ 
sulting  effective  interactions  are  obtained  by  averaging  over 
the  electronic  wavefunction  distribution  similar  to  the  case  of 
averaging  over  solvent  configurations.  Although  the  nuclear 
degrees  of  freedom  have  not  been  coarse-grained,  the  integra¬ 
tion  over  the  electronic  degrees  of  freedom  leaves  an  effective, 
“coarse-grained”  interaction  between  nuclei.  With  the  proper 
modification,  the  MS-CG  methodology  is  sufficiently  general 
and  provides  a  systematic  procedure  that  can  be  used  to  effi¬ 
ciently  parameterize  both  nonreactive  and  reactive  MD  mod¬ 
els,  such  as  the  framework  used  in  this  work,  that  account  for 
the  dynamic  formation  and  breaking  of  chemical  bonds. 

In  this  work,  a  new  iterative  FM  algorithm  is  developed 
that  enables  one  to  construct  reactive  models  using  flexible, 
tabulated  interactions  as  opposed  to  empirical  functions.  The 
resulting  models  enable  multiscale  reactive  molecular  dynam¬ 
ics  (MS-RMD)  models  for  chemical  reactions  in  condensed 
phase  systems  using  a  set  of  accurate  interaction  potentials. 
The  MS-RMD  method  has  the  additional  benefit  of  using  a 
linear  least-squares  method  to  determine  model  parameters, 
which  eliminates  the  complexities  associated  with  the  nonlin¬ 
ear  optimizations  of  empirical  functions. 

Section  II  will  first  summarize  the  AIMD  method  used  to 
generate  configurational  distributions  used  as  input  to  the  FM- 
based  MS-RMD  algorithm.  The  core  reactive  MD  methodol¬ 
ogy  is  then  briefly  reviewed.  Next,  attention  turns  to  nonre¬ 
active  models  for  both  the  hydronium  cation  and  hydroxide 
ion,  which  serve  as  the  foundation  for  the  corresponding  MS- 
RMD  models.  The  iterative  FM  algorithm  then  is  presented 
and  used  to  generate  reactive  models  for  both  charged-defects. 
Section  III  discusses  the  favorable  agreement  of  results  from 
the  FM  MS-RMD  models  with  those  from  the  original  AIMD 
simulations.  Section  IV  summarizes  the  work  and  its  potential 
benefits. 

II.  METHODS  AND  DEVELOPMENT 
A.  Condensed  phase  ab  initio  MD  simulations 

The  configurations  used  as  input  to  the  force-matching 
algorithms  discussed  in  this  work  were  generated  from 
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Born-Oppenheimer  molecular  dynamics  simulations  calcu¬ 
lated  by  using  the  Quickstep  module  in  the  CP2K  soft¬ 
ware  package.1  The  Becke-Lee- Yang-Parr  (BLYP)  (Refs.  32 
and  33)  density  functional  was  used  with  empirical  disper¬ 
sion  corrections34  for  all  AIMD  simulations.  The  use  of 
dispersion-corrected  density  functionals  lead  to  improved 
properties  of  water  at  ambient  conditions,  such  as  the  equilib¬ 
rium  density  and  an  increase  in  the  water  self-diffusion  con¬ 
stant  for  the  BLYP  functional.35-39  The  equations  of  motion 
for  the  nuclei  were  integrated  with  a  0.5  fs  timestep,  and  the 
electronic  orbitals  were  optimized  to  the  Born-Oppenheimer 
surface  by  using  an  orbital  transformation  method40  with  a 
convergence  criterion  of  1 0-7.  The  wavefunction  was  ex¬ 
panded  in  the  Gaussian  TZV2P  basis  set,  the  auxiliary  elec¬ 
tron  density  was  expanded  in  plane  waves  up  to  400  Ry,  and 
Goedecker-Teter-Hutter  pseudopotentials41  were  used.  The 
initial  configurations  for  both  the  excess  proton  and  hydroxide 
simulations  were  generated  by  adding  and  removing  a  proton, 
respectively,  from  a  previously  equilibrated  bulk  water  simu¬ 
lation  containing  128  water  molecules  with  a  lattice  constant 
of  15.66  A.  The  newly  created  systems  were  each  equilibrated 
in  the  constant  NVT  ensemble  for  25  ps  and  then  continued 
in  the  NVE  ensemble  for  an  additional  72  and  76  ps  for  the 
proton  and  hydroxide  systems,  respectively.  The  AIMD  sim¬ 
ulation  for  the  hydroxide  ion  discussed  here  is  the  same  one 
used  recently  to  develop  nonreactive  two- site  and  charged- 
ring  models  for  the  hydroxide  ion.21  In  addition,  two  extra 
32  ps  constant  NVE  simulations  were  completed  for  each 
charged  defect  in  order  to  calculate  dynamic  properties  and 
to  reduce  any  bias  of  results  due  to  initial  conditions.  The 
initial  configurations  for  these  additional  simulations  were  se¬ 
lected  from  simulations  calculated  with  the  corresponding  FM 
nonreactive  models  and  equilibrated  in  the  constant  NVT  en¬ 
semble  for  5  ps.  The  average  temperature  of  all  constant  ALE 
simulations  was  306.8  ±1.4  and  301.6  ±4.1  K  for  the  excess 
proton  and  hydroxide  systems,  respectively. 


B.  Reactive  MD 

The  modeling  of  chemical  reactions  requires  a  method¬ 
ology  that  can  dynamically  alter  the  bonding  topology  of 
the  system  over  the  course  of  a  simulation  in  response  to 
changes  in  the  environment  and  the  intrinsic  quantum  me¬ 
chanics  of  the  reaction.  Ab  initio  methods  that  explicitly  treat 
the  electronic  degrees  of  freedom  are  a  natural  choice  for  a 
methodology,  but  they  come  at  a  significant  computational 
expense.  One  set  of  solutions  to  address  this  challenge  that 
remain  feasible  even  for  large-scale  condensed  phase  systems 
are  multistate6,7,21'42  48  and  multiconfigurational49,50  molec¬ 
ular  mechanics  methods.  In  brief,  multistate  methods  treat 
the  system  as  a  linear  combination  of  several  possible  bond¬ 
ing  topologies  (diabatic  states)  that  are  coupled  to  one  an¬ 
other  through  the  off-diagonal  elements  of  a  quantum-like 
Hamiltonian  matrix.  In  particular,  the  MS-EVB  algorithm  is 
a  generalization  of  the  original  EVB  (Refs.  43  and  44)  ap¬ 
proach  that  adapts  and  dynamically  identifies  bonding  topolo¬ 
gies  to  include  as  the  simulation  progresses.  These  bonding 
topologies  form  a  basis  of  diabatic  states  that  are  used  to 
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FIG.  1.  Example  of  the  reactive  complexes  formed  for  the  (a)  excess  pro¬ 
ton  and  (b)  hydroxide  ion.  The  dashed  line  indicates  the  boundary  between 
those  atoms  within  the  reactive  complex  and  the  rest  of  the  system,  which 
dynamically  changes  over  the  course  of  a  simulation. 


evaluate  elements  of  the  Hamiltonian  matrix.  The  subset  of 
molecules  described  with  variable  bonding  topology  at  any 
point  in  time  during  a  simulation  form  the  reactive  complex. 
Examples  of  reactive  complexes  are  shown  in  Fig.  1  for  both 
the  excess  proton  and  hydroxide.  The  diagonal  entries,  Vjj , 
of  this  effective  Hamiltonian  typically  are  evaluated  by  using 
molecular  mechanics  forcefields,  Vmm,  with  additional  dia¬ 
batic  corrections,  Vcorr,  that  tend  to  be  repulsive.  The  off- 
diagonal  couplings,  Vjk,  provide  the  mechanism  for  the  sys¬ 
tem  to  transition  between  bonding  topologies,  and  therefore 
to  describe  chemical  reactions.  Typical  empirical  definitions 
for  the  off-diagonal  couplings  are  simple  geometric  functions, 
Ageo,  °f  the  donor-acceptor  distance,  such  as  Gaussian,  hy¬ 
perbolic  tangent  functions,  distributed  Gaussians,51  or  spline 
fits,18  to  reproduce  potential  energy  surfaces  calculated  from 
gas-phase  clusters.  Recently,  genetic  programs  have  also  been 
used  to  explore  both  function  and  parameter  space  for  defin¬ 
ing  interactions  in  MS-EVB  models.52  In  order  to  go  beyond 
just  the  local  geometry,  electrostatic  interactions,  Vex,  must 
be  included  in  the  condensed  phase  to  account  for  effects  due 
to  the  environments  that  may  alter  the  stability  of  configura¬ 
tions  near  a  transition  state.  The  specific  effective  Hamilto¬ 
nian  terms  are  given  by 

Vjj  —  Vmm  ±  Vcorr,  (1) 

Vjk  —  ( Vconst  ±  Vex)Ageo •  (2) 

When  electrostatic  interactions  are  included  in  the  off- 
diagonal,  a  constant,  Vconst,  that  would  normally  be  ab¬ 
sorbed  into  the  geometric  factor  as  a  prefactor  is  kept  separate 
in  order  to  shift  the  zero  of  the  electrostatic  energy.19- 42,47,48 
In  this  case,  the  electrostatic  energy  modulates  the  prefactor 
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of  the  original  geometric  factor.  The  diabatic  correction  term, 
Vcorr ,  used  here  was  labeled  in  previous  MS-EVB  models 
as  a  repulsive  interaction,  V/jep.19'24'42'47'48  This  term  is  la¬ 
beled  now  as  a  general  correction  because  no  constraints  on 
the  functional  form  are  enforced  in  the  FM  algorithm  pre¬ 
sented  here  and  it  is  thus  possible  that  this  interaction  could 
be  attractive  if  the  underlying  molecular  mechanics  forcefield 
had  too  large  of  an  energetic  penalty  for  sampling  regions  of 
phase  space  near  a  reactive  transition  state.  With  the  FM  algo¬ 
rithm  presented  here,  the  diabatic  corrections  and  off-diagonal 
geometric  factor  can  be  expressed  as  general,  flexible,  tabu¬ 
lated  potentials  as  opposed  to  a  set  of  empirical  functions  that 
may  or  may  not  have  been  validated.  As  such,  these  functions 
bridge  in  a  multiscale  fashion  the  MS-RMD  model  to  the  real 
quantum  Hellman-Feynman  forces  via  force-matching  into 
actual,  as  opposed  to  empirical,  effective  interactions.  In  the 
absence  of  empirical  functions,  the  interactions  in  the  result¬ 
ing  MS-RMD  models  can  be  used  to  justify  the  choice  of  em¬ 
pirical  functions  used  in  previous  work. 

C.  Force-matching  nonreactive  models 
1.  Hydronium  ion 

The  models  used  in  MS-RMD  simulations  are  based  in 
part  on  a  molecular  mechanics  forcefield,  which  serves  as  the 
major  contribution  to  the  diagonal  elements  of  the  Hamilto¬ 
nian  matrix  (we  note  that  this  is  not  in  principle  necessary  but 
it  is  convenient).  As  in  previous  work,  the  reactive  models  de¬ 
veloped  here  were  constructed  in  stages:  first  the  nonreactive 
model  was  developed,  and  then  the  reactive  model  was  con¬ 
structed.  The  Hellman-Feynman  theorem  (HFT)  forces  from 
the  cib  initio  Born-Oppenheimer  simulations  of  the  excess 
proton  were  used  as  input  for  a  FM  algorithm,  where  a  resid¬ 
ual,  Eq.  (3),  measuring  the  ensemble  average  of  the  squared 
difference  in  the  AIMD  HFT  forces  and  effective  forces  was 
minimized2 1 

x2=^(h*r-*r\2)-  (3> 

The  nonreactive  model  for  the  hydronium  cation  was 
parameterized  following  the  same  procedure  used  to  previ¬ 
ously  derive  a  model  from  a  condensed  phase  AIMD  simula¬ 
tion  using  the  Hamprecht-Cohen-Tozer-Handy  (HCTH)  den¬ 
sity  functional.24  The  intramolecular  model  for  the  hydronium 
cation  was  the  same  as  previously  used,  and  the  SPC/Fw  wa¬ 
ter  model53  was  used  as  opposed  to  fitting  a  new  water  model 
based  on  the  AIMD  simulation  results.  With  the  intramolecu¬ 
lar  and  electrostatic  parts  of  the  model  defined,  the  remaining 
nonbonded  interactions  between  the  hydronium  cation  and 
water  molecules  were  determined  by  force-matching  to  the 
HFT  forces  obtained  from  AIMD  simulations.  In  order  to  ac¬ 
complish  this,  a  subset  of  configurations  from  the  AIMD  tra¬ 
jectory  was  selected  for  the  FM  procedure,  such  that  the  solva¬ 
tion  structure  of  the  hydronium  cation  closely  resembled  the 
Eigen  cation,  HgO^,  the  resting  state  of  the  charged  excess 
proton  defect.54  55  The  criterion  for  selecting  these  configu¬ 
rations  was  based  on  the  smallest  value  of  a  proton  transfer 
order  parameter,  8  >  0.2  A,  calculated  as  the  difference  be¬ 


FIG.  2.  Pairwise  potentials  (solid  lines)  and  forces  (dashed  lines)  from  the 
B -spline  fits  for  hydronium- water  interactions  plotted  as  a  function  of  inter¬ 
atomic  separation.  The  atom  labels  are  defined  in  the  text. 


tween  the  H-bond  and  covalent  bond  distances  for  a  specific 
H-bonded  hydronium- water  pair.24,55 

The  interatomic  pairwise  potentials  to  be  fit  were  defined 
to  be  a  linear  combination  of  interpolating  functions  (cubic 
splines).  The  contributions  to  the  CG  forces  that  were  already 
defined  (intramolecular,  electrostatic,  SPC/Fw)  were  first  sub¬ 
tracted  from  the  HFT  forces,  and  the  atomic  force  residual 
was  used  in  the  minimization  of  Eq.  (3).  Each  interatomic 
pairwise  potential  was  represented  on  a  grid  with  a  spacing 
of  0. 1  A.  The  tabulated  force  and  its  first  two  derivatives  were 
constrained  to  be  zero  at  the  cutoff  distance  of  7.8  A.  The 
coefficients  for  these  spline  functions  were  determined  as  the 
linear  least-squares  minimization  of  the  residual  in  Eq.  (3)  by 
solving  an  overdetermined  system  of  equations.56-58  At  the 
shortest  distances  for  each  interatomic  pair,  where  the  sam¬ 
pling  is  limited,  the  forces  were  extrapolated  to  shorter  dis¬ 
tances  by  using  a  fit  of  1  If1,  where  the  exponent  n  was  chosen 
based  on  smoothness.  These  extrapolated  force  curves  were 
then  fit  to  B-splines  in  order  to  eliminate  statistical  noise  in  the 
force  curves  and  improve  energy  conservation.  The  B-spline 
fitted  force  curves  for  the  hydronium-water  interactions  are 
shown  in  Fig.  2,  and  the  coefficients  are  listed  in  Table  I.  The 
atom  types  used  in  Fig.  2  and  elsewhere  are  labeled  as  two 
characters,  with  the  first  letter  defining  the  element  and  the 
second  letter  defining  the  molecule.  For  example,  “OH”  de¬ 
notes  the  oxygen  of  hydronium  (or  hydroxide  depending  on 
the  model)  and  “HW”  corresponds  to  the  hydrogen  atom  of  a 
water  molecule. 

2.  Hydroxide  ion 

The  nonreactive  hydroxide  model  developed  here  is 
based  on  the  two-site  model  from  a  previous  study21  with 
an  additional  three-body  interaction  included.  Hydroxide  ions 
commonly  are  over-coordinated  in  bulk  water  when  using 
simple  two-site  models,  and  thus  more  complicated  mod¬ 
els  are  required  using,  for  example,  three-body  potentials  or 
a  charge-ring  model.1516'20'21  With  the  two-site  hydroxide 
model21  as  a  starting  point,  a  three-body  interaction  involv¬ 
ing  the  hydroxide  bond  and  proximal  water  oxygen  atoms 
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TABLE  I.  B-spline  coefficients  for  each  hydronium-water  interatomic  pair 
in  the  nonreactive  model.  The  values  of  rcore  used  in  the  B-spline  fits  are  2.0, 
1.2,  0.9,  and  1.2  A  for  the  OH-OW,  OH-HW,  HH-OW,  and  HH-HW  inter¬ 
atomic  pairs,  respectively.  The  evenly  spaced  bin  widths  for  the  correspond¬ 
ing  curves  were  0.58,  0.33,  0.23,  and  0.33  A.  The  value  of  rcut  remained  at 
7.8  A  and  a  polynomial  order  of  6  was  used  for  all  curves.  Trailing  zeroed 
coefficients  for  each  curve  are  omitted  for  clarity,  and  the  resulting  curves  are 
shown  in  Fig.  2. 


n 

OH-OW 

OH-HW 

HH-OW 

HH— HW 

i 

142.3787 

30.6048 

158.1527 

13.1760 

2 

119.7021 

12.8241 

76.6490 

3.7248 

3 

74.9776 

-  1.3372 

16.3827 

-3.5881 

4 

45.3455 

-  7.5645 

1.4867 

-3.4542 

5 

-  20.9843 

-  5.9280 

-2.2173 

-  6.0079 

6 

11.2329 

-  9.0470 

-2.7144 

-4.3451 

7 

-3.9884 

3.6121 

-3.5742 

0.1435 

8 

1.5298 

-4.5015 

-5.2106 

0.2671 

9 

-3.3958 

2.1488 

5.0724 

-  1.1518 

10 

1.0618 

-  6.2932 

-  4.0822 

-0.8201 

11 

-  1.1024 

3.1139 

4.5130 

-  1.0234 

12 

0.9083 

-0.8365 

1.5914 

-0.3816 

13 

-  0.4423 

1.2911 

3.8124 

-  1.8704 

14 

0.8274 

0.4376 

-0.1383 

-0.0519 

15 

-0.3076 

0.3791 

1.5763 

-  1.4096 

16 

0.1095 

0.0871 

1.5115 

-  0.3603 

17 

-0.0013 

2.1562 

-  0.5254 

18 

0.7607 

1.5052 

-0.2139 

19 

-0.3421 

1.8256 

-0.3881 

20 

1.2020 

1.4905 

-  0.4447 

21 

-  0.0946 

1.2609 

-  0.3272 

22 

0.1670 

1.1820 

-  0.0429 

23 

-0.2593 

1.0116 

-0.4616 

24 

0.1893 

0.6997 

-0.1153 

25 

-0.0451 

1.0547 

-  0.0973 

26 

0.0228 

0.5954 

0.0153 

27 

0.4978 

28 

0.7561 

29 

0.6394 

30 

0.2425 

31 

0.3365 

32 

0.3206 

33 

0.3029 

34 

0.1039 

35 

0.0483 

36 

-  0.0025 

(HH-OH-OW)  was  included  in  the  model.  The  optimal  pa¬ 
rameters  for  a  Stillinger- Weber  (SW)  potential,59  Eq.  (4), 
were  determined  by  using  a  Simplex  algorithm  and  fitting  the 
residual  force  difference,  Eq.  (3),  defined  between  the  original 
two-site  hydroxide  model  and  the  A1MD  simulation  results. 
The  form  of  the  SW  potential  used  is  given  by 

Vsw(n,  r2,0)  =  s(cosd— cosdo)2  exp  ( - )exp( - )• 

\rx-a  J  \r2-o  J 

(4) 


The  optimized  values  of  the  parameters  are  s  =  271.6 
kcal/mol,  6  =  -0.290  rad,  and  a  =  4.0  A.  As  discussed  be¬ 
low,  the  addition  of  this  three-body  potential  to  the  nonreac¬ 
tive  model  produces  a  considerable  improvement  in  the  hy¬ 


droxide  ion  solvation  structure,  as  observed  in  the  hydroxide- 
water  angle  distribution  (HH-OH-OW). 


D.  Force-matching  reactive  models 
1.  Iterative  FM  algorithm 


A  straightforward  application  of  a  linear  least-squares 
method  for  the  parameterization  of  a  MS-RMD  model  is  com¬ 
plicated  by  the  fact  that  the  forces  are  known  only  after  the 
multistate  Hamiltonian  has  been  constructed  and  diagonalized 
and  the  Hellman-Feynman  theorem  used  with  the  coefficients 
of  the  ground-state  eigenvector,  such  that 


3H 


=  -('t'ctl  —  |4/°)  = 


'inn 

j  ’ 


(5) 


If  the  coefficients  for  the  MS-RMD  ground  state  for 
an  ensemble  of  configurations  were  known  a  priori,  then  it 
would  be  straightforward  to  express  each  component  of  the 
force,  F )mn,  as  a  linear  combination  of  interpolating  func¬ 
tions  (cubic  splines),  and  the  resulting  total  atomic  force,  Fy, 
would  be  a  linear  function  of  the  model  parameters  to  be  de¬ 
termined.  However,  this  is  not  the  case.  Another  possible  is¬ 
sue  is  the  nonlinear  dependencies  introduced  in  the  definition 
of  the  MS-RMD  model,  such  as  the  off-diagonal  couplings 
being  defined  as  an  electrostatic  energy  multiplied  by  a  lo¬ 
cal  geometric  factor,  Eq.  (2).  If  the  geometric  factor,  Aqeo, 
was  defined  to  be  a  linear  function  of  model  parameters,  then 
the  nonlinearity  introduced  by  multiplication  with  the  electro¬ 
static  energy  and  constant  could  be  handled  by  using  two  sep¬ 
arate  optimizations:  linear  and  nonlinear.  With  an  initial  guess 
for  the  nonlinear  parameters,  a  linear  least-squares  optimiza¬ 
tion  could  be  used  to  determine  those  parts  of  the  model  that 
depend  linearly  on  model  parameters,  such  as  Aqeo •  After  the 
linear  optimization  completes,  the  nonlinear  parameters  could 
be  optimized  holding  the  new  set  of  linear  parameters  fixed. 
Afterwards,  the  linear  portion  of  the  model  could  be  reopti¬ 
mized  by  using  the  new  set  of  nonlinear  parameters  and  the 
procedure  repeated  until  converged.  By  defining  a  large  frac¬ 
tion  of  the  model  as  being  linear  with  respect  to  the  model 
parameters,  the  cost  of  the  nonlinear  optimization  can  be  sig¬ 
nificantly  reduced  by  requiring  the  optimization  of  only  a  few 
parameters. 

If  one  chooses  to  use  this  type  of  iterative  algorithm  to 
reduce  the  cost  of  the  parameter  optimization  due  to  nonlin¬ 
earities  in  the  model,  then  an  immediate  solution  to  the  prob¬ 
lem  of  not  knowing  the  ground  state  MS-RMD  coefficients  a 
priori  is  to  provide  an  initial  guess  and  then  include  the  con¬ 
vergence  of  the  eigenvector  as  part  of  the  criteria  for  reaching 
self-consistency.  This  is  the  main  motivation  behind  the  FM 
algorithm  outlined  in  Fig.  3  for  a  single  set  of  values  for  the 
nonlinear  parameters.  In  this  algorithm,  an  initial  Hamiltonian 
matrix  is  constructed  by  using  educated  guesses  for  all  model 
parameters,  the  Hamiltonian  matrix  is  diagonalized,  and  the 
coefficients  of  the  ground  state  eigenvector  determined.  With 
this  set  of  eigenvector  coefficients  and  fixed  values  for  the 
nonlinear  parameters,  the  atomic  forces  in  Eq.  (5),  which  are 
linear  functions  of  the  remaining  parameters,  are  variationally 
optimized  by  using  a  linear  least-squares  method  to  minimize 
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FIG.  3.  Flow  chart  outlining  the  iterative  force-matching  algorithm  used  to 
develop  reactive  models. 


the  residual  in  Eq.  (3).  With  the  new  set  of  linear  parameters,  a 
new  guess  for  the  Hamiltonian  matrix  can  be  constructed  and 
diagonalized  in  order  to  update  the  coefficients  of  the  ground 
state  eigenvector  to  be  used  to  calculate  a  new  set  of  linear 
parameters.  This  procedure  is  iterated  until  self-consistency 
has  been  reached  for  the  ground  state  eigenvector  and  linear 
model  parameters.  This  algorithm  for  defining  the  MS-RMD 
interactions,  Vcorr  and  Ageo,  as  linear  combination  of  in¬ 
terpolating  functions  was  found  to  converge  in  three  to  five 
iterations  for  both  the  proton  and  hydroxide  ion  models  using 
forces  calculated  from  an  AIMD  condensed  phase  simulation. 
If  there  were  nonlinear  parameters  in  the  model  to  be  deter¬ 
mined,  they  would  be  optimized  by  using  an  algorithm  such  as 
Simplex  with  the  other  portions  of  the  model  fixed,  including 
the  linear  parameters,  identical  to  the  procedure  used  in  previ¬ 
ous  work.24  With  the  current  algorithm,  however,  the  number 
of  nonlinear  parameters  has  been  significantly  reduced.  The 
procedure  described  here  for  separately  optimizing  the  linear 
and  nonlinear  parameters  can  then  be  repeated  sequentially 
until  both  sets  of  parameters  have  converged  to  their  optimal 
values. 

The  main  advantages  of  this  algorithm  are  the  ability  to 
define  interactions  as  flexible  tabulated  expressions  and  the 
reduction  in  computational  cost  associated  with  a  nonlinear 
optimization.  For  the  reactive  models  developed  here,  the  re¬ 
maining  nonlinear  parameters  are  contained  in  the  definition 
of  the  off-diagonal  coupling  as  the  prefactor  to  the  geometric 
factor.  The  off-diagonal  constant  and  exchange  charges  are 
the  only  four  parameters  that  need  to  be  determined  by  using 
a  nonlinear  optimization  algorithm.  The  diabatic  corrections 
and  local  geometric  factor  are  described  by  using  spline  func¬ 
tions,  the  same  definition  used  for  the  pairwise  interactions  of 
the  nonreactive  hydronium  and  hydroxide  models. 

2.  Excess  proton 

For  the  reactive  model  describing  the  hydrated  proton, 
the  initial  guess  for  the  four  nonlinear  parameters  was  taken 
from  a  previously  developed  model  derived  from  condensed 
phase  simulations  with  the  HCTH  density  functional.24  Sim¬ 


ilar  to  the  MS-EVB3  model  for  the  excess  proton,  two  di¬ 
abatic  corrections  were  defined  as  part  of  the  model  devel¬ 
oped  here:  one  for  the  hydronium-oxygen  and  water-oxygen 
atoms,  VPPP,  and  one  between  hydronium-hydrogen  and 
water-oxygen  atoms,  VP°P,  such  that 

Vcorr  —  Vrep  +  VP°P.  (6) 

The  off-diagonal  geometric  factor,  Ageo ,  was  defined  to 
be  a  function  of  the  distance  between  the  hydronium-oxygen 
and  water-oxygen  involved  in  the  H-bond  of  the  transferring 
proton.  The  forces  for  each  of  these  three  pairwise  func¬ 
tions  were  represented  on  a  grid  with  a  spacing  of  0.1  A. 
The  tabulated  force  and  its  first  two  derivatives  were  con¬ 
strained  to  be  zero  at  the  cutoff  distance  of  2.9,  1.3,  and  3.2 
A  for  the  Vggp,  VPPP,  and  Ageo  interactions,  respectively. 
For  the  initial  iteration,  the  diabatic  corrections  were  set  to 
zero,  and  the  off-diagonal  geometric  factor  was  set  to  one 
for  all  donor-acceptor  distances.  With  the  iterative  algorithm 
outlined  in  Fig.  3,  the  coefficients  for  these  spline  functions 
were  determined  as  the  linear  least-squares  optimization  of 
the  residual  in  Eq.  (3)  by  solving  an  overdetermined  system 
of  equations.56-58  At  the  shortest  distances  for  each  pairwise 
function,  the  forces  were  extrapolated  to  shorter  distances  by 
using  a  linear  fit.  These  extrapolated  force  curves  were  then 
fit  to  B-splines  to  eliminate  statistical  noise  in  the  force  curves 
and  improve  energy  conservation.  The  B-spline  fitted  curves 
were  used  for  each  step  of  the  iterative  FM  algorithm  in  Fig.  3. 
The  B-spline  fitted  force  curves  for  the  reactive  proton  model 
are  shown  in  Fig.  4,  and  the  polynomial  coefficients  are  listed 
in  Table  II.  In  the  final  set  of  four  nonlinear  parameters,  the 
off-diagonal  exchange  charges  were  found  to  be  similar  to  the 
initial  values  used  at  the  start  of  the  fitting  process.  Based 
on  this  result,  the  iterative  FM  algorithm  was  repeated  with 
the  off-diagonal  exchange  charges  fixed  at  the  values  from 


r(A) 


FIG.  4.  Tabulated  potentials  (solid  lines)  and  forces  (dashed)  for  diabatic 
corrections  (a)  VPPP,  (b)  Vppp  in  Eq.  (6)  and  off-diagonal  geometric  factor 
(c)  Ageo  in  Eq.  (2)  for  the  excess  proton  (black)  and  hydroxide  ion  (red) 
reactive  models. 
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TABLE  II.  B-spline  coefficients  for  the  diabatic  corrections  and  off- 
diagonal  geometric  factor  for  the  reactive  excess  proton  model.  The  values 
of  rcore  used  in  the  B-spline  fits  are  2.0,  0.8,  and  2.1  A  for  the  Vppp,  VPPP, 
and  Aqeo  interactions,  respectively.  The  evenly  spaced  bin  widths  for  the 
corresponding  curves  were  0.085,  0.045,  and  0.0525  A.  The  value  of  rcu,  was 
2.85,  1 .25,  and  3. 15  A  for  the  corresponding  curves,  and  a  polynomial  order 
of  6  was  used  for  all  curves.  Trailing  zeroed  coefficients  for  each  curve  are 
omitted  for  clarity,  and  the  resulting  curves  are  shown  in  Fig.  4. 


n 

yOO 
y  REP 

yHO 
v  REP 

AqEO 

i 

186.8016 

41.2160 

2.7165 

2 

183.1398 

40.3683 

2.7125 

3 

175.7715 

38.6700 

2.6192 

4 

164.7942 

36.1281 

2.6283 

5 

150.0428 

32.7296 

2.4218 

6 

131.7579 

28.4960 

2.4187 

7 

109.6329 

23.3915 

2.2148 

8 

87.7021 

18.3660 

1.8722 

9 

65.5203 

13.3516 

3.3914 

10 

43.7042 

8.7389 

2.9196 

11 

23.6742 

4.6277 

3.0227 

12 

10.3514 

2.1909 

2.9306 

13 

3.8574 

0.5109 

2.8816 

14 

1.0625 

0.3833 

2.4445 

15 

0.4695 

0.0700 

2.3205 

16 

0.1798 

0.0486 

2.6821 

17 

3.1176 

18 

3.1232 

19 

2.8403 

20 

2.0784 

21 

1.2542 

22 

0.5621 

23 

0.2077 

24 

0.0695 

25 

0.0269 

26 

0.0069 

the  previous  model,  leaving  only  the  off-diagonal  constant  as 

a  parameter  to  be  optimized  by  using  a  simple  line  search. 

The  values  for  the  nonlinear  parameters  for  the  reactive  pro- 

ton  model  are  listed  in  Table  III. 

3.  Hydroxide  ion 

The  development  of  a  reactive  model  for  the  hydrox- 

ide  ion  followed  the  same  procedure  as  discussed  in  Sec.  II 

D  2  for  the  excess  proton.  The  definition  of  this  model  in- 

eluded  similar  definitions  for  the  diabatic  corrections  and  off- 

diagonal  geometric  factor  as  pairwise  functions  between  the 

hydroxide-oxygen  and  water  molecule  atoms.  The  initial  val- 

TABLE  III.  Optimized  model  parameters  for  the  off-diagonal  couplings  for 

FM  reactive  proton  and  hydroxide  models. 

Excess  proton 

Hydroxide 

-  13.88 

-  13.78 

«£(e) 

-  0.0500 

-0.0510 

q“S  (e) 

0.0167 

0.0247 

9&W 

0.0332 

0.0526 

ues  for  the  four  nonlinear  parameters  were  again  taken  from 
the  previous  excess  proton  model,24  as  opposed  to  selecting 
a  random  set  of  initial  values.  The  details  of  fitting  the  three 
tabulated  potentials  to  spline  functions  were  the  same  as  for 
the  reactive  proton  model  developed  in  Sec.  II  C.  Although 
the  use  of  an  HH-OH-OW  three-body  correction  was  found 
to  improve  the  solvation  structure  of  the  nonreactive  two-site 
hydroxide  model,  initial  tests  with  reactive  models  developed 
with  this  correction  were  found  to  yield  long  hydroxide  life¬ 
times  and  thus  a  low  rate  of  successful  proton  transfer  events. 
Motivated  by  the  three-body  interaction  used  in  a  recently  de¬ 
veloped  reactive  model  for  the  hydroxide  ion,15  a  three-body 
interaction  was  defined  to  act  between  the  hydroxide  oxy¬ 
gen  and  two  neighboring  water  oxygens  (OW-OH-OW)  us¬ 
ing  as  a  basis  the  same  parameters  as  for  the  nonreactive  HH- 
OH-OW  interaction.  While  the  OW-OH-OW  interaction  was 
not  sufficient  to  reduce  the  degree  of  over-coordination  ob¬ 
served  in  the  underlying  nonreactive  two-site  model,  initial 
tests  with  reactive  models  were  promising  and  thus  this  was 
the  interaction  form  used  to  develop  the  model  presented 
here.  If  one  wanted  to  choose  a  more  flexible  basis,  then 
additional  three-body  corrections  could  be  force-matched 
using  AIMD  data,  such  as  the  coordination-dependent  (3, 
4,  and  5  solvating  water  molecules)  OW-OH-OW  interac¬ 
tion  used  in  a  recent  MS-EVB  model  for  the  hydroxide 
ion.15 

Results  from  an  initial  reactive  hydroxide  model  showed 
an  increased  population  of  the  5-coordinated  solvation  struc¬ 
ture,  with  a  peak  in  the  hydroxide-water  angle  (HH-OH- 
OW)  probability  distribution  near  180°,  which  was  origi¬ 
nally  observed  in  the  nonreactive  model  before  introduction 
of  the  three-body  potential.  With  the  current  definition  for 
the  local  geometric  factor,  which  depends  only  on  the  donor- 
acceptor  distance,  there  is  no  distinguishing  between  water 
molecules  in  the  square  planar  arrangement  around  the  hy¬ 
droxide  ion  or  a  water  molecule  that  approached  the  hydrox¬ 
ide  ion  from  below.  Since  transfer  of  a  proton  from  a  wa¬ 
ter  molecule  that  approaches  below  the  hydroxide  ion  was 
not  observed  in  the  AIMD  simulations,  the  off-diagonal  ge¬ 
ometric  factor  was  multiplied  with  an  additional  function 
that  depends  on  the  HH-OH-OW  angle.  This  function  was 
chosen  to  reduce  the  off-diagonal  coupling  when  a  water 
molecule  donates  an  H-bond  at  an  HH-OH-OW  angle  near 
180°,  while  leaving  the  coupling  for  all  other  water  molecules 
unchanged, 

fm=  ^(l-tanh[«(0-0„)]).  (7) 

The  values  of  the  parameters  were  chosen  to  be  a 
=  0.25°- 1  and  00  =  140°  based  on  the  HH-OH-OW  angle 
distribution  as  calculated  from  the  AIMD  simulations,  Fig.  5. 
This  neglect  of  certain  diabatic  states  could  have  also  been 
implemented  as  part  of  the  state  search  algorithm,  where  an 
angle  criterion  would  be  used  to  determine  which  states  are 
included  in  the  diabatic  basis  used  to  construct  the  Hamilto¬ 
nian  matrix,  but  that  would  have  likely  resulted  in  states  with 
significant  coupling  abruptly  being  included  or  removed  from 
the  Hamiltonian,  leading  to  an  increase  in  the  drift  of  the  total 


Downloaded  15  Aug  2012  to  130.202.50.2.  Redistribution  subject  to  AIP  license  or  copyright;  see  http://jcp.aip.org/about/rights_and_permissions 


22A525-8  Knight,  Lindberg,  and  Voth 


J.  Chem.  Phys.  137,  22A525  (2012) 


FIG.  5.  Probability  distributions  for  the  hydroxide  water-oxygen  angle  (HH- 
OH-OW)  for  the  AIMD  simulation  (solid  lines)  and  FM-based  MS-RMD 
models  (dashed  lines).  The  full  distribution  (black)  is  decomposed  into  con¬ 
tributions  from  solvation  structures  with  three  (red),  four  (green),  and  five 
(blue)  water  molecules  donating  H-bonds  to  the  hydroxide  ion,  with  relative 
populations  indicated. 


energy.  By  using  a  smooth  function  of  the  system  coordinates, 
the  contribution  from  those  states  where  a  water  molecule  is 
donating  an  H-bond  from  below  the  hydroxide  ion  can  be 
made  negligible  in  a  continuous  manner.  With  the  parameters 
for  the  angle-dependent  function  fixed,  the  iterative  FM  pro¬ 
cedure  was  repeated  in  order  to  parameterize  a  reactive  model 
for  the  hydroxide  ion.  The  optimized  nonlinear  parameters  are 
listed  in  Table  III,  the  coefficients  for  the  tabulated  potentials 
are  listed  in  Table  IV,  and  the  tabulated  potentials  are  shown 
in  Fig.  4. 


TABLE  IV.  B-spline  coefficients  for  the  diabatic  corrections  and  off- 
diagonal  geometric  factor  for  the  reactive  hydroxide  model.  The  values  of 
r core  used  in  the  B-spline  fits  are  2.0,  0.9,  and  2.3  A  for  VPPp .  Vg°p,  and 
Aqeo  interactions,  respectively.  The  evenly  spaced  bin  widths  for  the  corre¬ 
sponding  curves  were  0.085,  0.035,  and  0.085  A.  The  value  of  rcut  was  2.85, 
1.25,  and  3.15  A  for  the  corresponding  curves  and  a  polynomial  order  of  6 
was  used  for  all  curves.  Trailing  zeroed  coefficients  for  each  curve  are  omit¬ 
ted  for  clarity  and  the  resulting  curves  are  shown  in  Fig.  4. 


n 

V°0 

v  REP 

yHO 
v  REP 

Ageo 

i 

333.9791 

6.1891 

5.9091 

2 

327.4629 

6.0462 

5.9784 

3 

314.5365 

5.7586 

5.6280 

4 

294.9678 

5.3307 

5.9403 

5 

269.1527 

4.7566 

5.2448 

6 

236.5075 

4.0525 

5.2644 

7 

197.7903 

3.2104 

4.8854 

8 

158.5645 

2.4098 

5.7094 

9 

120.0611 

1.6668 

3.8874 

10 

80.2545 

1.0334 

3.7608 

11 

43.4640 

0.5137 

4.0839 

12 

18.4555 

0.2308 

1.7564 

13 

6.8357 

0.0921 

0.9765 

14 

1.8423 

0.0046 

0.0154 

15 

0.7601 

-  0.0074 

0.1343 

16 

0.2668 

-  0.0229 

-0.0106 

III.  RESULTS  AND  DISCUSSION 

A.  Force-matched  reactive  models 

The  tabulated  potentials  for  the  diabatic  corrections  and 
off-diagonal  geometric  factor  are  shown  in  Fig.  4  and  found 
to  take  similar  shapes  for  both  reactive  proton  and  hydrox¬ 
ide  models.  The  diabatic  corrections  take  a  relatively  sim¬ 
ple  shape  and  are  repulsive  for  all  distances  sampled.  This 
result  is  in  agreement  with  previous  MS-EVB  models  that 
have  defined  these  corrections  as  repulsive  interactions  mod¬ 
eled  by  using  an  exponential  function.16,  ly. 42,47,48  por  Gff_ 
diagonal  geometric  factors,  however,  there  is  structure  in  the 
forces,  F qeo,  as  a  function  of  the  donor-acceptor  distance  for 
an  H-bond  that  the  proton  is  transferring  across,  but  the  cor¬ 
responding  geometric  factors  are  relatively  featureless  when 
viewed  at  the  scale  in  Fig.  4.  These  peaks  in  the  derivatives  of 
the  geometric  factor  would  be  missed  with  simple  Gaussian  or 
hyperbolic  tangent  functions  typically  used  in  EVB  models. 
The  functional  form  for  the  geometric  factor  used  in  the  MS- 
EVB  models  for  the  excess  proton,19, 42,47,48  which  requires 
the  optimization  of  seven  parameters,  would  also  have  diffi¬ 
culty  capturing  the  structure  observed  in  the  range  of  sampled 
distances  2.25-3.2  A. 

The  force-matched  reactive  models  developed  in  this 
work  were  used  to  generate  independent  trajectories  in  the 
constant  NVE  ensemble  by  using  a  modified  version  of  the 
LAMMPS  MD  code.60  The  equations  of  motion  for  all  sim¬ 
ulations  with  FM  models  (nonreactive  and  reactive)  were 
integrated  by  using  a  timestep  of  0.5  fs  unless  otherwise 
stated.  The  state  search  algorithm  developed  for  the  MS- 
EVB3  model  was  used  here  with  generalizations  to  model  the 
transport  of  both  the  excess  proton  and  hydroxide  ion.19 

The  total  energy  drift  for  the  nonreactive  FM  models  was 
0.16  ±  0.24  and  0.26  ±  0.13  kcal/mol  per  nanosecond  for 
the  hydronium  and  hydroxide  models,  respectively,  as  esti¬ 
mated  from  six  2.5  ns  simulations  each.  The  energy  drifts  rel¬ 
ative  to  the  total  energy  were  9.2  x  10~5  and  1.3  x  1 0  4 
1/ns  for  the  nonreactive  hydronium  and  hydroxide  models. 
The  total  energy  drift  for  the  reactive  proton  model  was  3.3 
±  3.8  kcal/mol  per  ns  (2.1  x  10-3  1/ns),  as  estimated  from 
sixteen  400  ps  simulations  using  a  timestep  of  0.50  fs.  The 
drift  in  the  total  energy  for  the  hydroxide  reactive  model  was 
5.8  ±  4.8  kcal/mol  per  ns  (3.4  x  10~3  1/ns),  as  estimated 
from  twenty-two  300  ps  simulations  using  a  timestep  of  0.5 
fs.  With  a  smaller  timestep  of  0.25  fs,  the  drift  for  the  hy¬ 
droxide  model  was  reduced  to  3.0  ±  3.8  kcal/mol  per  ns  (1.8 
x  10  3  1/ns).  The  results  of  simulations  for  both  reactive 
models  are  discussed  below. 

B.  Excess  proton 

The  hydronium-water  radial  distribution  functions 
(RDFs)  for  the  AIMD,  nonreactive,  and  reactive  models  are 
shown  in  Fig.  6.  Comparison  of  the  four  RDFs  shows  a  gen¬ 
eral  agreement  with  the  AIMD  data  for  both  models  in  terms 
of  the  peak  heights,  positions,  and  integrated  coordination 
numbers.  One  noticeable  deviation  is  in  the  RDF  for  the 
hydronium-oxygen  water-oxygen  (OH-OW)  interatomic  pair 
with  an  enhanced  peak  at  separations  just  larger  than  3.0  A. 
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r(A) 


FIG.  6.  Radial  distribution  functions  (solid  lines)  and  integrated  coordination 
numbers  (dashed  lines)  between  atoms  of  the  hydronium  cation  and  water 
molecules.  For  each  plot,  results  from  the  AIMD  (black),  nonreactive  (blue), 
and  reactive  MS-RMD  (red)  models  are  shown.  The  atom  labels  are  defined 
in  the  text. 


The  density  contributing  to  this  peak  principally  comes  from 
the  second  solvation  shell  peak  (near  4.5  A),  where  a  water 
molecule  that  transitions  in  and  out  of  the  first  solvation 
shell  and  which  would  be  expected  to  weakly  donate  an 
H-bond  to  the  lone  pair  side  of  the  hydronium  cation  is 
bound  too  strongly.  As  discussed  below  for  the  hydroxide 
ion,  it  seems  to  be  a  general  trend  in  these  reactive  models  for 
charged  defects  that  special  attention  needs  to  be  focused  on 
correctly  describing  relatively  weak  interactions  compared 
with  the  strong  H-bonds  formed  in  the  first  solvation  shell. 
The  inclusion  of  a  second  set  of  pairwise  potentials  that  are 
zero  beyond  the  first  solvation  shell  may  assist  in  accurately 
capturing  features  in  the  solvation  structure  due  to  weaker 
H-bonds.  There  are  some  deviations  in  the  RDFs  beyond 
the  second  solvation  shell,  but  this  can  be  expected  because 
the  forces  in  that  region  are  dominated  by  water-water 
interactions  and  the  SPC/Fw  water  model  used  differs  from 
the  corresponding  AIMD  water  model. 

An  important  thermodynamic  property  for  the  transport 
of  charged  defects  via  the  Grotthuss  mechanism  is  the  free 
energy  barrier  for  proton  transport.  The  proton  transfer  barrier 
calculated  from  the  AIMD  simulation  is  0.72  kcal/mol  for  the 
CP2K  simulation,  Fig.  7,  which  is  larger  than  the  barriers  of 
~  0.6  and  0.5  kcal/mol  calculated  for  the  BLYP  and  HCTH 
density  functionals,  respectively.24,61  The  barrier  height  for 
the  FM  reactive  model  is  0.08  kcal/mol  smaller  than  the 
AIMD  result.  The  increased  barrier  height  for  proton  trans¬ 
port  with  respect  to  the  other  density  functionals  is  consis¬ 
tent  with  a  corresponding  reduction  in  the  diffusion  constant 
of  the  excess  proton.  The  diffusion  constant  for  the  ex¬ 
cess  proton  calculated  from  the  current  AIMD  simulations  is 
0.44  ±  0.33  A2/ps  to  be  compared  with  the  diffusion  constants 
of  ~  0.07  and  0.31  A2/ps  for  the  BLYP  (Ref.  62)  and  HCTH 
(Ref.  24)  functionals.  Diffusion  constants  of  0.07  ±  0.06  and 
0.45  ±  0.20  A2/ps  were  calculated,  respectively,  for  the  non¬ 
reactive  and  reactive  FM  models  developed  here,  where  the 
value  for  the  nonreactive  model  corresponds  to  solely  vehicu¬ 


FIG.  7.  Potentials  of  mean  force  along  a  proton  transfer  coordinate  for  the 
excess  proton  and  hydroxide  models.  Results  are  shown  for  the  AIMD  (solid 
black),  nonreactive  (dashed  blue),  and  reactive  (dot-dashed  red)  models. 


lar  transport  (classical  diffusion).  The  experimental  value  for 
the  proton  diffusion  constant  is  0.94  A2/ps.63  These  calculated 
values  for  the  proton  diffusion  coefficients  are  considerably 
smaller  than  the  experimental  value,  but  this  can  be  explained 
by  the  slower  self-diffusion  of  the  water  combined  with  the 
fact  that  these  are  classical  simulations  and  nuclear  motion 
has  not  been  quantized.'11  The  self-diffusion  of  water  for  the 
current  AIMD  simulation  is  0.04  ±  0.01  A2/ps,  which  is  in 
agreement  with  a  recently  reported  value  using  the  same  func¬ 
tional  and  dispersion  correction. 39  The  self-diffusion  of  water 
in  the  SPC/Fw  water  model  used  to  develop  the  FM  models 
for  the  excess  proton  is  0.23  ±  0.05  A2/ps,  which  agrees  with 
the  experimental  value.53,64 

C.  Hydroxide  ion 

For  the  nonreactive  FM  hydroxide  model,  the  inclusion 
of  the  three-body  Stillinger- Weber  potential  leads  to  improved 
agreement  with  the  AIMD  results  compared  with  results  of 
the  original  two-site  model.21  The  agreement  with  the  cur¬ 
rent  model  in  describing  the  hydroxide  ion  solvation  structure, 
shown  in  Fig.  5,  is  similar  to  that  of  the  charged-ring  model, 
which  was  also  force-matched  to  AIMD  data,21  as  determined 
by  comparison  of  the  HH-OH-OW  angle  distributions.  The 
corresponding  reactive  model  also  agrees  with  the  AIMD  re¬ 
sults,  such  that  the  main  peak  in  the  full  distribution  is  near 
107°.  There  is  still  a  noticeable  contribution  close  to  180°  in 
the  reactive  FM  model,  which  largely  arises  from  4-  and  5- 
coordinated  species.  The  population  of  5-coordinated  species 
for  the  nonreactive  model  is  larger  than  the  AIMD  result  and 
the  population  is  further  enhanced  in  the  reactive  FM  model 
with  a  corresponding  decrease  in  the  3-coordinated  species. 
The  relative  fraction  of  4-coordinated  species  in  the  nonre¬ 
active  and  reactive  FM  models  is  nearly  identical  and  ~12% 
smaller  than  the  AIMD  result. 

The  RDFs  for  both  the  nonreactive  and  reactive  hydrox¬ 
ide  models  are  shown  in  Fig.  8.  The  same  feature  observed 
just  beyond  the  first  peak  in  the  OH-OW  RDF  for  the  hydro¬ 
nium  models  can  also  be  seen  in  the  RDFs  for  the  hydroxide 
models.  Again,  there  is  an  enhanced  peak  between  the  two 
main  peaks  observed  in  the  AIMD  RDF,  resulting  from  water 
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FIG.  8.  Radial  distribution  functions  (solid  lines)  and  integrated  coordina¬ 
tion  numbers  (dashed  lines)  between  atoms  of  the  hydroxide  ion  and  water 
molecules.  For  each  plot,  results  from  the  AIMD  (black),  nonreactive  (blue), 
and  reactive  MS-RMD  (red)  models  are  shown.  The  atom  labels  are  defined 
in  the  text. 


molecules  being  held  too  tightly.  In  this  case,  it  is  due  to  the 
water  molecules  that  act  as  H-bond  acceptors.  The  increased 
binding  is  consistent  with  increased  peak  heights  and  inte¬ 
grated  coordination  numbers  for  the  two  hydroxide-oxygen 
RDFs.  Additional  evidence  for  the  increased  binding  of  water 
molecules  that  accept  an  H-bond  from  the  hydroxide  ion  can 
be  seen  in  the  HH-OW  RDF  and  the  narrowness  of  the  first 
peak  near  2.0  A  for  the  reactive  model.  The  corresponding 
peak  for  water  molecules  that  accept  an  H-bond  from  the  hy¬ 
droxide  ion  is  not  present  in  the  nonreactive  model,  which  is 
likely  a  result  of  the  three-body  HH-OH-OW  potential  acting 
on  all  nearby  water  molecules. 

The  free  energy  barrier  for  proton  transfer  between  water 
and  hydroxide  is  1 .6  and  1 .9  kcal/mol  for  the  AIMD  and  reac¬ 
tive  models,  respectively,  as  seen  in  Fig.  7.  The  self-diffusion 
of  the  hydroxide  ion  in  the  AIMD  simulation  was  calculated 
to  be  0.05  ±  0.02  A2/ps,  which  is  larger  than  the  calculated 
water  self-diffusion  constant,  0.02  ±  0.01  A2/ps,  and  smaller 
than  that  of  the  excess  proton  discussed  above.  This  ordering 
of  the  diffusion  constants  is  consistent  with  previous  AIMD 
results  and  experiment.65  The  hydroxide  diffusion  constants 
calculated  from  the  FM  models  are  0.08  ±  0.05  and  0.28 
±  0.20  A2/ps  for  the  nonreactive  and  reactive  models,  re¬ 
spectively.  The  diffusion  constant  for  the  reactive  MS -RDF 
model  developed  here  is  similar  to  the  value  calculated  from 
a  charged-ring  MS-EVB  model,  0.31  ±  0.03  A2/ps.16 


IV.  SUMMARY  AND  FUTURE  OUTLOOK 

In  this  work,  a  new  force-matching-based  reactive  MD 
method  was  discussed  with  the  goal  of  no  longer  relying  on 
predefined  empirical  functions  for  the  reactive  models.  With 
this  iterative  algorithm,  numerical  tabulated  potentials  can  be 
used  for  all  interactions  including  diabatic  corrections  and 
off-diagonal  couplings  in  the  MS-RMD  framework.  The  flex¬ 
ibility  and  generality  of  the  algorithm  are  illustrated  by  the 
development  of  reactive  models  for  the  excess  proton  and 


hydroxide  in  bulk  water,  which  are  nontrivial  systems  to 
model.  The  FM  procedures  used  in  this  work  are  independent 
of  the  source  of  the  reference  data,  thus  enabling  one  to  incor¬ 
porate  data  from  a  number  of  configurational  sampling  meth¬ 
ods  and  increasingly  higher  levels  of  theory.  Using  a  general 
procedure  to  develop  reactive  models  with  a  reduced  number 
of  constraints  on  the  definition  of  the  model,  one  can  gen¬ 
erate  increasingly  more  accurate  models  that  can  reproduce 
the  physical  properties  of  the  reference  AIMD  (or  related) 
method,  while  at  the  same  time  extending  the  time  and  length 
scales  accessible  to  the  simulations  at  a  significantly  reduced 
computational  cost. 
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